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Some fundamental problems for an energy conserving adaptive resolution molecular 

dynamics scheme 

Luigi Delle Sitf] 

Max- Planck- Institut fur Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany 

Adaptive resolution molecular dynamics (MD) schemes allow for changing the number of degrees of 
freedom on the fly and preserve the free exchange of particles between regions of different resolution. 
There are two main alternatives on how to design the algorithm to switch resolution using auxiliary 
"switching" functions; force based and potential energy based approach. In this work we show that, 
in the framework of classical MD, the latter presents fundamental conceptual problems which make 
unlikely, if not impossible, the derivation of a robust algorithm based on the potential energy. 
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Introduction 



Multiscale modeling and simulation in condensed mat- 
ter is a field of continuous expansion as the basic prop- 
erties of an increasing number of systems, relevant to 
current research, are discovered to strongly depend on a 
delicate scales' interplay. The massive progress of com- 
puter technology together with the parallel development 
of novel powerful simulation methods has strongly con- 
tributed to this expansion so that by now detailed se- 
quential studies from the electronic scale to the meso- 
scopic and continuum are routinely performed. However 
among all these methods of a particular interest are those 
which deal in a more direct way with the multiscale idea. 
Typically these schemes are based on a single computa- 
tional approach which links two or more interconnected 
scales. One example is the technique used to study edge 
dislocation in metals, where local chemistry affects large 
scale material properties or the crack of materials 
where the rupture of a local interatomic bond is then 
propagated to the larger scale and again back to the next 
interatomic bond and so on 0, Q ; m this case quantum 
based methods are interfaced with classical atomistic and 
continuum models within a single computational scheme. 
A further example is the Quantum Mechanics/Molecular 
Mechanics scheme [H . This is mainly used for soft matter 
systems where a fixed region of space requires quantum 
resolution and the external part is treated at classical 
atomistic level. Examples are solvation of large molecules 
where the chemistry happens locally (quantum region) 
while the statistical effect of the fluctuating environment 
(solvent) far from the molecules can be treated in a rather 
efficient way at classical level. In the same fashion there 
are several more examples (see e.g. Refs.fH @). All of 
these anyway are characterized by a non trivial limita- 
tion, i.e. the region of resolution is fixed and free ex- 
change of particles with the other regions are not allowed. 
While this may not be a crucial point for system involv- 
ing rigid structures, certainly is a very strong limitation 
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for highly fluctuating systems. The natural next step 
to overcome this problem is the design of novel adap- 
tive resolution methods which indeed allow for the ex- 
change of particles among regions of different resolution. 
In general, in such a scheme a molecule moving from a 
high resolution region to a lower one, would gradually 
loose some degrees of freedom until the lower resolution 
is reached but yet the statistical equilibrium among the 
two different regions is kept at any instant. Recently 
some schemes based on this idea have been presented 
in the literature @, 1, They differ in the way the 
different resolutions are coupled in the MD algorithm. 
The coupling can be achieved either through the poten- 
tial, slowly passing from an atomistic to a corresponding 
coarse grained potential (and vice versa), or through the 
forces, that is slowly passing from a force derived from 
an atomistic potential to a force derived from the corre- 
sponding coarse grained potential (and vice versa). The 
passage from the atomistic to the coarse grained is con- 
trolled by a smooth "switching function" which is used 
to interpolate the two quantities. For the force based 
scheme it is not possible to define a potential energy from 
the interpolation formula, but on the basis of physical ar- 
guments this problems can be circumvented (lot [Tl| as it 
will be briefly discussed later on; for the potential based 
scheme obviously the definition of potential energy is the 
central point. In this sense, the potential based scheme 
would seem more appealing, however the subject of this 
work is to show that on the basis of a mathematically rig- 
orous derivation, this scheme is not applicable. Here we 
construct the most general adaptive scheme based on the 
potential and derive the necessary conditions by which 
one can obtain the switching functions. As an outcome 
we show that the resulting set of partial differential equa- 
tions has got boundary conditions such that the system 
is overdetermined and thus solutions may exists only for 
trivial cases. Moreover, even in case a solution may exist, 
further technical problems, due to the nature of the dif- 
ferential equation, arise which make this scheme rather 
unpractical. The paper is organized as follows; in the 
next section a short overview of the force based method 
is presented, it summarizes its crucial point and enumer- 
ates the latest applications. Next the potential based 
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scheme is presented with its general features. Finally, a 
general interpolation scheme is used to derive the equa- 
tion that defines the switching functions. The paper is 
closed by the discussion and conclusions. 

B. Force based scheme: A short overview of the 
AdResS method 

According to the previous discussion, a method which 
has turned to be rather robust is the Adaptive Resolution 
Simulation (AdResS) @, [1] ■ It is based on coupling the 
atomistic and the mesoscale through an interpolation for- 
mula for the atomistic and coarse grained force. At this 
point it should be mentioned that this approach, as well 
as the calculations performed in this work, are valid, so 
far, under the assumption of pair interactions; however 
due to the large use of pair potentials in atomistic simu- 
lation, the AdResS method, as well as the result of this 
work are nevertheless of interest to the simulation com- 
munity. Briefly, the space is divided in two regions as for 
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FIG. 1: (Color on line) Schematic picture of the partitioning 
of space in high resolution (atomistic) region B, low resolution 
(coarse grained) region A and transition region A. w(x) is 
the switching function which allows a smooth transition from 
a coarse grained to an atomistic resolution and vice versa. 
Below the pictorial representation of a tetrahedron molecule 
that change resolution according to the position in space is 
presented. This representation is taken from [3] 

example in FigQ] a high resolution region, let us call it B 
where the molecule has atomistic resolution and a region 
A where the molecule is coarse grained. In between there 
is a region A where a smooth transition from one resolu- 
tion to another takes place via a continuous "switching" 
function w(x), such that w(xi) = 0;w(x 2 ) — 1- The 
interpolation formula then reads 0] : 

Fq/j = w(X a )w(X[3)F^p m + [1 — w(X a )w(Xp)\F™p (1) 

where a and (3 labels two distinct molecules, F^ m is 
derived by the atomistic potential where each atom of 
molecule a interacts with each atom of molecule /3, and 
F^g is obtained from an effective pair potential between 
the centers of mass of the coarse-grained molecules; the 
latter is derived on the basis of the reference all-atom 



system. Eq[I] does not allow to define a potential in the 
switching region [ToL fllj , however, in this scheme such 
a definition is not required. Actually, all one needs to 
know is based on the following arguments: the change of 
resolution can be interpreted in terms of similarity with 
a geometrically induced first order phase transition with 
an associated latent heat This interpretation justi- 
fies the use of a thermostat, during an MD simulation, 
in the switching region A so that the physical equilib- 
rium is kept. Numerical calculations and applications to 
rather different systems have shown that indeed this ap- 
proach is satisfactory (see the applications to a liquid of 
tetrahedral molecules Q, a polymer solvated in it, [ijj 
and liquid water [III). The crucial point of this scheme 
is that EqfT] is only an ansatz based on satisfying the 
third Newton law and on numerical simplicity; however 
the numerical results show that the method indeed gives 
the correct answers when compared with all atom simu- 
lations and its physical interpretation is consistent with 
the basic principles of equilibrium in statistical mechan- 
ics [lH . However, one may naturally ask whether the 
same or a similar interpolation scheme can be applied to 
potentials and thus preserve the energy conservation as 
suggested by Ensing et al. @. In the following section 
we show that to build an interpolation scheme similar to 
EqfT] but applied to potentials instead offerees it is not 
possible. 



C. A generic scheme based on the potentials 

Let us define two generic switching functions: 
f(X a ,Xp) continuous and differentiablc in A, and out- 
side A defined such that: 

f(X a ,Xp) — 0; X a > x 2 ,and,Xf3 > x 2 

f(X a ,Xp) = l;X a ,or,Xp< Xl (2) 

and g(X a , Xp) continuous and differentiable in A and 
outside A defined such that: 

g{X a , Xp) — ; X a < xiorXp < x\ 

g(X a , Xp) = 1; X a andXp > x 2 . (3) 

Here X a and Xp are the coordinates along the x direc- 
tion, as represented in Fig[TJ of the center of mass respec- 
tively of the generic molecule a and f3. A generalization 
of EqfT] to the potentials using these two generic switch- 
ing function f{x) and g{x) writes: 

jjcoupUng = f {Xa ,Xp)U cg + 9 {X a , Xp)U atom (4) 

where U coup u ng is the potential coupling the two resolu- 
tions, U cg = t/ cg (R Q ,R^) is the coarse grained poten- 
tial and R a ,Hp the coordinates of the centers of mass; 
Uatom = U a tom{ Y aii r 0j) ls the atomistic potential be- 
tween atom i of molecule a and atom j of molecule 
P. EqQ] couples the different scales similarly to what 
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is done by EqQ] but the with the hypothetical advan- 
tage of automatically conserving energy. At this point 
to do molecular dynamics, we need to derive the forces. 
The following situations are clear: if the molecules are 
located both in region A (coarse grained force), or both 
in region B (atomistic force), or one in A and one in B 
(coarse grained force), or one in A and one in A (coarse 
grained force). However once the molecules are both in 
A or one in A and one in B, the force must be de- 
rived by the whole expression of Eq0] Let us calculate 
the coupling force acting on R Q and R^. One should 
keep in mind that R Q = J2i=i n r ai/ n an d equivalently 
Hp = , / ri, where for simplicity the molecules 

where chosen to have both n identical atoms. It follows 
that the force acting on the center of mass of molecule a 
is : 



aTJcoupling 

^coupling 



<9R Q 



(5) 



which in explicit form writes: 



^upUng = _ f{x ^ Xp) ^E. g(Xa ,Xp)^^ - 

8f(X a ,X f3 ) dgjX^Xp) 

Uc.a U atom ^ 1,0 J 



'eg 

taking into account that 
as: 



dR 



<9R Q 



^y 2 - = 1, Eq[6]can be rewritten 



F coupling ft v 
r„ - n A o 



X fi )F r 



■g(X a ,Xp)F a c 



■ drift 



where 



(7) 



drift 



-u. 



df(X a ,Xt 
' dX a 



dg{X ai X p ) 

Uatom KTP I s J 

OA a 



is a spurious force with no physical meaning which 
emerges as a consequence of the presence of the switch- 
ing functions. In fact in this case the center of mass 
of a molecule receives an additional acceleration in the 
switching region, due to the switching function, which 
should not be there because from the physical point of 
view the molecules in any resolution regime must be 
equivalent and the switching function is not a physical 
quantity. The effect that it will have is of drifting par- 
ticles along the x direction. At this point, the condition 
to have a well based physical treatment of the particles 
without artifacts in the dynamics due to the introduction 
of the switching functions is to recover from Eq[71 in a 
mathematical way, the force coupling scheme involving 
only the atomistic and coarse grained force. This will al- 
low us to determine f(X a ,Xp) and g(X a , Xp) for which 
the energy is conserved and to have an algorithm that 
in principle works rather well as shown by the AdResS 
scheme. The physical condition to do so, as implicitly 
suggested by Ensing et ai, is: 



drift 







(9) 



and in this case, translated into the mathematical condi- 
tion, becomes: 



U c 



df{X ai Xp) 
' dX a 



U at 



dg(X a ,Xp) 
1 dX a 



Q 



(10) 



To make the problem mathematically correct one should 
follow the same procedure for the force acting on R^g so 
that the final conditions reads: 



u c 



df(x a ,Xf)) 

9 dX a 
df(X a ,Xp) 



dg{X a ,Xp) n 

U atom ^ , r U 



dX, 



U at 



dX a 
dg[X a ,X p ) 







dX 



0. 



(11) 



This is a system of first order partial differential equa- 
tions where g and / are the unknown functions in A and 
Xfj are the variables 



14j. Without going into the 



X a , 

details of the mathematical properties of such a system, 
a simple and yet powerful observation clearly shows that 
a solution may exist only in very special cases but cer- 
tainly not in general. This observation is rather simple; a 
differential equation or a system of differential equations 
of the first order has got solutions which are uniquely 
identified by one boundary condition (one for / and one 
for g in this case). At this point if one goes back to 
the definition of / and g given in Eqs l2l31 it is easy to 
see that in order to have a valid switching function with 
the correct limiting case at the boundary of A, there are 
two boundary conditions for each function, associated to 
Eafm to be satisfied: 



f(X a ,Xp) = 0; X a = x 2 ,and,Xp = x 2 
f(X a ,Xp) — 1; X a — xi,and,Xp = x\ 



(12) 



and 



g(X a , Xp) — 0; X a = xiand, Xp = x\ 
g(X a ,Xp) = 1; X a = x 2 ,and,Xp = x 2 - (13) 

This means the system of equations is overdetermined 
and a solution in general does not exist . Specifically, 
if the equation are solved using the condition in xi, it 
may or may not exist a solution such that / in a cer- 
tain point x 2 is equal to zero, and equivalently for g; 
of course the same arguments is valid if as a boundary 
condition is chosen that of x 2 . However even in case a 
solution exists, there would be no control on the switch- 
ing region A as it is not possible to locate one of the two 
boundaries a priori. This aspects makes this approach 
not convenient for any robust MD algorithm. A further 
point that invalidates the potential approach is the fact 
that, while ideally / and g should be function solely of 
X, to deal with a simple algorithm, Eg I 111 shows that in- 
deed at least one of the two functions should depend on 
all the degrees of freedom of the atomistic system, as in 
the equation the atomistic potential depends on all such 
degrees of freedom. This may be even possible for simple 
systems, however as the molecules become larger this ap- 
proach becomes highly unpractical. One may even think 
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of a more general scheme in the same fashion of what 
is proposed in Ref.[§], that is to introduce an additional 
potential $ such that the coupling potential reads: 

jjcoupHng = f ^ + g ^ X )U atom + $ (14) 

where $ is equal to zero in A and B , in order to ob- 
tain U cg in A and U a tom in B and it is a certain regular 
function in A such that: 



' drift 



<9$ 

dX~ 



a, (3; V X a ,Xp e A. 



(15) 



In this way one obtains a more general expression for the 
energy in the switching region and, regarding the forces, 
the role of Jy- is that of removing the spurious force 
due to the switching functions. At this point one must 
notice that EqlTS] is equivalent to Eq[TU] (or Eq fTT|) with 
the only difference of the presence of Jy- on the r.h.s. 
The conclusions of this work do not change because the 
problem of the boundary conditions of this differential 
equation remains the same as for Eq[TU] (or Eq fTT|) . In 
this case however there is more flexibility and one can dis- 
tinguish two situations: (a) $ is a known function and 
g and / unknown; (b) $ is unknown and g and / are 
known. In case (a) the conclusions drawn before do not 
change because we still have two first order partial dif- 
ferential equations in g and / and the overdetermination 
is not removed by the presence of the known term — -j^- 
on the r.h.s of Eq[ll] In case (b) we will have again a 
system of first order partial differential equations where 
the unknown function is $ (Eq [To| and is characterized 
by two boundary conditions, one in x\ and one in X2, 
(i.e. $ = 0), thus the overdetermination is shifted from 
/ and g to 



D. Conclusions 



We have shown that an adaptive resolution method 
based on the ansatz of potential interpolation via switch- 
ing functions cannot be realized as the mathematical con- 
dition of finding a suitable switching function it is likely 
to have no solution or only trivial ones. It was already 
shown before that this scheme leads to the violation of 
Newton third law [l(J El for the special case / = 1 — g. 
The arguments presented above add up to the previous 
one and further show that for the most generic interpola- 
tion formula the switching functions do not exist except 
for some special and trivial cases. In general, for a nu- 
merical implementation, such a scheme would not be fea- 
sible. This fact does not exclude the possibility that the 
adaptive resolution can be achieved via other approaches 
based on the potential energy. In fact, recently Hyden et 
al. 15j have presented an alternative scheme for adap- 



tive resolution based on potentials. This scheme, rather 
promising, can be applied also to the quantum-classical 
interface. However it does not make use of switching 
functions and looses the numerical simplicity of the in- 
terpolation formula together with its physical interpreta- 
tion which instead is the non trivial advantage of the force 
based method. In conclusion, the development of adap- 
tive resolution approaches is a field of rapidly growing 
interest, the intention of this work is that of fixing some 
clear directions along which one can or cannot move in 
order to develop more sophisticated and yet numerically 
simple schemes. 

I am grateful to Matej Praprotnik and Kurt Kremer 
for helpful comments on the manuscript. This work is 
supported in part by the Volkswagen foundation. 
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